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Abstract 


The vapor flow in a heat pipe was mathematically modeled and the 
equations governing the transient behavior of the core were solved numerically. 
The modeled vapor flow is transient, axisymmetric (or two dimensional) 
compressible viscous flow in a closed chamber. The two methods of solution 
are described. The more promising method failed (a mixed Galerkin finite dif- 
ference method) whereas a more common finite difference method was suc- 
cessful. Preliminary results are presented showing that multi-dimensional 
flows need to be treated. 

A model of the liquid phase of a high temperature heat pipe was developed. 
The model is intended to be coupled to a vapor phase model for the complete 
solution of the heat pipe problem. The mathematical equations are formulated 
consistent with physical processes while allowing a computationally efficient 
solution. The model simulates time dependent characteristics of concern to the 
liquid phase including input phase change, output heat fluxes, liquid temper- 
atures, container temperatures, liquid velocities and liquid pressures. Prelimi- 
nary results were obtained for two heat pipe startup cases. The heat pipe 
studied used lithium as the working fluid and an annular wick configuration. 
Recommendations for implementation based on the results obtained are pre- 
sented. 

Experimental studies were initiated using a rectangular heat pipe. Both twin 
beam laser holography and laser Doppler anemometry were investigated. Pre- 
liminary experiments were completed and results are reported. 


xii 



Chapter I 

INTRODUCTION 

A heat pipe is a device which transfers heat by evaporating and condensing a 
fluid. A wick structure is used to return the liquid to the evaporation zone by 
capillary effect. Figure 1 shows a typical sketch of a heat pipe which consists 
of the case, the wick structure and the vapor space. The input heat at one end 
of the pipe increases the temperature of the liquid in the wick structure. The 
liquid then evaporates, the pressure in the vapor core increases and the vapor 
flows to the other end. At this end the temperature is lower and the vapor 
condenses and releases heat. The liquid, then, flows from the condenser to the 
evaporator by capillary action through the wick structure. 

The heat pipe is considered as a heat transfer device with very high 
conductivity. It can operate in weightless environments and transfer a huge 
amount of heat under conditions of very low temperature differences. These 
advantages make the heat pipe unique in space applications, electronic cooling 
systems and fusion energy systems. Although the steady state operation of 
heat pipe is well established, the transient behavior is not entirely well under- 
stood. 

The reliability of a heat pipe for space applications depends on its ability 
to respond to a sudden change in the heat load, and startup and shut down 
of the thermal system. A comprehensive experimental and theoretical study is 
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Figure 1: A Heat Pipe 

needed to understand and properly model the physical phenomena which oc- 
cur in heat pipes during operational transients. 

The objectives of this work were to investigate the dynamic behavior of the 
vapor flow in heat pipes, experimentally and theoretically. Experiments were 
to be carried out using a planar heat pipe to examine the transient response 
of the vapor core to a sudden increase in the input heat flux. Optical methods, 
which do not interfere the physical process, were proposed for measurements 
of velocity and temperature Fields. The experimental apparatus and the optical 
methods are described in Chapter V. 

The vapor flow in a heat pipe was mathematically modeled and the 
equations governing the transient behavior of the vapor core were solved nu- 


2 



mcricallv. The modeled vapor flow is a transient, axisymmctric- (or two- 
dimensional). compressible, viscous flow in a closed chamber. The boundary 
conditions include inflow and outflow of heat and mass which simulate the 
boundaries of the vapor core in a heat pipe. The proposed numerical method, 
along with preliminary results, are discussed in Chapter VI. The liquid phase 
is considered in Chapter VII. 

The proposed objectives of the study were not fully met. Computer models 
were developed only to the state where feasibility of an approach was demon- 
strated. An experimental apparatus was built and put through some prelimi- 
nary shake down testing. Much remains to be done to satisfy the objectives. 
This, of necessity, is left to others. 
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Chapter II 

APPLICATIONS OF HEAT PIPES 

Gaugler [1] was the first to develop the concept of a heat pipe in 1942 pro- 
posing to apply heat pipes for cooling the interior of an ice box. Later in 1964, 
an intensive work in this field was started by Grover and his coworkers [2] 
at the Los Alamos National Laboratory. The objective was the application 
of the heat pipe in spacecraft designs. The first quantitative analysis of heat 
pipes was performed by Cotter in 1965 [3]. Since then there have been many 
advances in theory, design and practice of heat pipes. 

The fact that heat pipes are able to operate under very low gravity condi- 
tions and that they do not include any moving parts, has motivated the devel- 
opment of the heat pipe for space applications. In these applications the 
temperature is very high (1000- 2000“C ) and liquid metals (such as sodium, 
lithium) are used as working fluids. 

A heat pipe was proposed by Grover [2] to supply heat to the emitters of 
thermionic electrical generators and to remove heat from the collectors of these 
devices. In space based nuclear power systems, heat pipes are considered as a 
means to transfer heat from the generating point to the radiators and also to 
transfer heat throughout the radiators, providing a uniform temperature on 
the radiator surfaces [4], 
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Recently heat pipes have been proposed to be used in internal cooling of the 
leading edges of reentry vehicles and hypersonic aircrafts to reduce the peak 
temperature and alleviate the thermal gradients at the edge [5]. In this case 
the working fluid (a liquid metal) which is initially in the solid state, is melted 
by the reentry heat load and the heat pipe operation is initiated. The use of 
high temperature heat pipes as a heat transfer device has also been considered 
in rocket engines to cool and make the nozzles isothermal [6]. 

In case of a loss of coolant accident in a fusion reactor, it is desirable to 
ensure that no first wall melting or evaporation occurs. Such melting or evap- 
oration of constituent elements can lead to the dispersal of radioactive waste 
hazards. In contrast to fission systems, a fusion power core usually contains 
enough thermal capacity to distribute the after-heat efficiently. Since the ma- 
jor component of decay after-heat is within a thin zone close to the first wall 
structure, a design incorporating heat pipes would result in an effective and 
fast redistribution of the decay after-heat. A heat pipe with a proper working 
fluid may be chosen such that the mechanism of heat transport through the 
pipe will start above a threshold design temperature. This would make the 
heat pipe operative only when needed, that is in the case of a loss of coolant 


accident. 



Chapter III 

THEORETICAL BACKGROUND 

Theory of heat pipe operation is briefly outlined in this section. Pressure drop 
and heat flux in a heat pipe are discussed for engineering calculations along 
with the limits on the heat pipe operation. The section is ended with discussion 
of the transient operations of heat pipes. The detailed explanation of heat pipe 
theory may be found in Dunn and Reay [7], 


3.1 PRESSURE DROP 

In order for heat pipes to operate, the capillary driving force in the wick, 
A p cap , must overcome the pressure drop in the liquid phase, A p t , in the vapor 
phase, A p v , and the gravitational head, if any, A p g . That is 

bPcap > Api + A p v + A p g ( 1 ) 

The pressure difference across the interface of liquid and vapor is 


where a is the surface tension, 6 is the contact angle, and r is the radius of 
curvature. 

Figure 2 illustrates the operational conditions of the vapor-liquid interface. 
The radius of curvature in the condenser section is very large, whereas evapo- 
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ration causes the radius of curvature at the evaporator mcniscii to decrease. 
The difference between the radii of curv ature in condenser and evaporator, 
provides the capillary force. Thus, the capillary driving force, using Eq. 2, is 


^Pcap 




( 3 ) 


where the subscripts cand c refer to evaportator and condenser, respectively, 
and r is the radius of the capillary pore. Usually r c is very large with respect 
to r e . This yields 


^ 'Pcap 


2 a cos 6 


(4) 


Vapor 



Figure 2: Operational Level of Interface Between Liquid and Vapor 

Pressure drop in the liquid phase, A p t , depends mainly on the structure of 
the wick and the properties of the liquid. For engineering calculations the 
liquid phase is modeled as a liquid flow in a porous medium. Since the mass 
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flow in the evaporator and condenser is not constant, an effective length. l eff . 
is defined for the porous medium. The effective length, through which the 
mass flow is constant, is approximated as 

l e ff= l a + 

where l a is the length of the adiabatic section which is between the evaporator 
and condenser sections. l e and l c are the lengths of the evaporator and 
condenser. Thus, Darcy's law is used to relate the pressure drop in the liquid 
phase to the wick and liquid properties, i.e., 


Api = 


Pi K 


( 6 ) 


where and p t are viscosity and density of the liquid, m is the liquid flow rate, 
A w is the wick cross sectional area and K is the wick permability given by 


K — 


2er h 

fl Re l 


where e is the wick porosity, r h is the hydraulic radius and f, Re t is the friction 
factor-Reynolds number product for laminar duct flow. The values of K for 
different types of wick structure is tabulated in Chi [8]. 

The pressure drop in the vapor phase, A p v , consists of the pressure drop in 
the adiabatic section, A p va , and the pressure drop in the evaporator, A p ve and 
in the condenser, A p vc . That is, 

A p Y = A p va + A p ve + A p vc (7) 
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In order to calculate for the pressure drop in the adiabatic section, the vapor 
flow in this region is assumed to be laminar Poiseuille flow in a circular tube 
[7], Then the pressure drop A p Ya , is given by 


A Pva 


8p v m l a 

4 

n p v r v 


(8) 


where r v is the radius of the vapor core. Applying a one-dimensional analysis 
Cotter [3] calculated the pressure drop in the evaporator and condenser, 
which are formulated as 


A Pve = 


m 


ve a 

8 P v r* 


a - m 

A Pvc = 3 4 

2w“p v r v 


Upon substituting Eqs. (8) and (9) in Eq. (7), we get 


8 p v mL m 2 

A Pv = ~ ~ ~4~ 0.074 

P v r v p v r v 


(9 a) 


(9 b) 


( 10 ) 


However, Busse [9], by using a modified Poiseuille velocity profile in Navier- 
Stokes equations, has shown that A p v can be expressed as 


A p v = 


8p v m l e ff 

4 

^ Pv 


( 11 ) 


where l eff is defined in Eq. (5). 
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The hydrostatic head. A p a . is written as 

&p g = g Pi L Sin (p (12) 

where <p is the inclination angle of the heat pipe and L is the total length of the 
heat pipe. 

Upon substituting Eqs. (6), (11) and (12) into Eq. (1), the circulation rate, 
m , is given by 


kPcap ~ gP/ L sin (p 


m = 


V/ 


Pl KA w 


+ 


8/* v 

4 

KPSv 


(13) 


3.2 HEAT FLUX 

Neglecting the sensible heat which is small compared to the latent heat, the 
axial heat flux rate, q, is given by 

q = rhhf g (14) 

where, m is the circulation rate of the working fluid and hj s is the latent heat 
of evaporation. 

From Eqs. (13) and (14) the heat transfer rate can be expressed as a func- 
tion of the capillary driving force, the geometry and the properties of the 
working fluid; i.e., 
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P = 


(APcap ~ SPl L sin <P) h 


l fg 


e fJ\ 


Pi 


Pl KA w 


+ 


8 n v 

4 

np v r v 


(15) 


3.3 LIMITATIONS 

There are several limitations on the output heat transported through heat 
pipes. These limitations are briefly described below. 

1. Viscous limit: In long pipes and at low temperatures, the vapor pressure is 
low and the effect of viscous friction on the vapor flow may dominates over the 
inertial forces. In this situation the circulation of the working fluid is limited, 
which consequently, limits the heat transfer through the pipe. 

2. Sonic limit: At low vapor pressures, the velocity of the vapor at the exit of 
the evaporator may reach the speed of sound. Then the evaporator cannot 
respond to further decrease in the condenser pressure. That is, the vapor flow 
is chocked, which limits the vapor flow rate. 

3. Capillary limit: A capillary structure is able to provide circulation of a given 
fluid up to a certain limit. This limit depends on the permability of the wick 
structure and the properties of the working fluid. 

4. Entrainment limit: The vapor flow exerts a shear force on the liquid in the 
wick which flows opposite the direction of the vapor flow. If the shear force 
exceeds the resistive surface tension of the liquid, the vapor flow entrains small 
liquid droplets (Kelvin-Hclmholtz instabilities). The entrainment of liquid in- 
creases the fluid circulation but not the heat transfer through the pipe. If the 
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capillary force cannot accommodate the increased flow, dryout of the wick in 


evaporator may occur. 

5. Boiling limit: At high temperatures, nucleate boiling may take place which 
produces vapor bubbles in the liquid layer. The bubbles may block the wick 


pores and decrease the vapor flow. Furthermor: 


presence of bubbles de- 


creases the conduction of heat through the liquid layer which limits the heat 


transfer from the heat pipe shell to the liquid which is by conduction only. 



Figure 3: Heat Pipe Operating Limits 
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The above limits are fully explained in several textbooks (for example, Chi 
[8]) and are not reviewed here. These limitations dictate the operational en- 
velope for a given design of a heat pipe, as shown diagramaticallv in Fig. 3 
[2]. The shape of the operational envelope changes for different working flu- 
ids and different wick materials. 


3.4 TRANSIENT OPERATION 

Transient operation of heat pipes may be considered in terms of three different 
regimes: startup, operating transients and shut down. The following is a brief 
description of these regimes. 

Sta rtup 

Heat pipe startup is a transient operation through which the heat pipe starts 
its operation from an initial static condition. The working fluid of heat pipe 
may initially be very cold or frozen in the wick. Under this condition, the va- 
por pressure is very low and the vapor flow is free molecular. The heat 
transfered from the evaporator melts the frozen working fluid in the 
evaporator and evaporation takes place in the liquid-vapor interface. The va- 
por pressure is initially very low and the vapor flow is highly accelerated. Then 
the flow may be chocked at the exit from the evaporation zone if its velocity 
reaches the sonic velocity. By addition of heat, the frozen working fluid is 
completely melted and a continuum vapor flow prevails throughout the vapor 
core. 
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Behavior of the vapor flow during the startup of a heat pipe from a frozen 
state, may be considered in three distinct phases: (/) free molecular flow 
throughout the vapor core, (ii) establishment of continuum flow in the evapo- 
ration zone with a front, moving towards the condenser with possibility of 
chocked flow at the evaporator exit and (Hi) continuum flow of vapor all 
through the vapor core. 

Ogerating Transients 

Operating transients involve power changes during nominal operating con- 
ditions. There arc several transient conditions of interest. A heat pipe may 
operate in a low power stand-by mode but be required to quickly respond to 
reach full power operation. A heat pipe may also be required to perform under 
load-following conditions, which depend on changing demands of the power 
conversion system. A nuclear reactor control system may also require changes 
in heat generation that would have to be handled by the heat pipe system. 

The physics of interest for these situations depend on the initial conditions, 
the rate of change and magnitude of input power and condenser conditions. 
A heat pipe operating in stand-by mode at low vapor pressure may encounter 
the sonic limit while powering up. Large temperature gradients and complex 
vapor dynamics would develop. After passing through the sonic-limited oper- 
ation, or for a heat pipe initially operating at high vapor pressure, the 
capillary, entrainment and boiling limits are of concern. 

Shut Down 

Shut-down of a heat pipe refers to the transient process occurring when the 
power of the thermal system is shut-off. The input heat of the heat pipe is 
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suddenly reduced to zero and the pipe cools down to a static condition. There 
are several problems encountered with this mode of transient operation. After 
the power is shut off. the working fluid of the heat pipe cools down and even- 
tually solidifies. This phase change may affect the wick structure. The screen 
sheets of the wick are very delicate structures with fine meshes. Freezing of the 
working fluid in the screen holes may damage the wick structure by tearing the 
screen. Also if for some reason, the working fluid freezes mainly inside the 
screen holes, it will cause a major reduction in the driving force of the liquid 
flow. 

Another possible problem is the appearance of cracks or bubbles in the 
solidified mass of the working fluid. The presence of cracks or bubbles can 
decrease the conduction of heat in the solidified mass during the startup, 
thereby the subsequent startup time increases. Cracks can also cause strain 
on the structure of the heat pipe and the screen and fatigue may become a se- 
rious problem after many shut downs. A proper heat pipe design must take 
into account the problems involved in the shut down transient mode, in order 
to have a reliable and long lasting heat pipe. 



Chapter IV 

LITERATURE REVIEW 


The heat pipe studies were initiated by the experimental work of Grover et al. 
[2] and the quantitative analysis of Cotter [3]. Initial qualitative experiments 
were performed on heat pipes with water and sodium as working fluids. Cotter 
considered a right circular heat pipe with large length-to-diameter ratio. He 
performed a simple one-dimensional analysis on the vapor flow to develop the 
general basic theory for quantitative engineering calculations of heat pipe op- 
erations. These studies have been continued to develop the theory and per- 
formance of heat pipes. 

The vapor flow in heat pipes is a complicated problem due to high nonlinear 
character of the governing equations and the inflow and outflow boundaries 
in the evaporator and condenser. Different approaches have been used to 
simplify the problem. In most of the work done so far, the vapor flow is ana- 
lyzed under steady state conditions, as a one dimensional flow [10,11,12,13] 
or a two dimensional flow [14-24]. 

Banskton and Smith [14,15] developed numerical solutions for 
incompressible, steady state, two dimensional vapor flow in a heat pipe. Their 
results show flow reversal in the condenser for radial Reynolds numbers larger 
than two. Indication of a reverse flow is important in design of a heat pipe in 
calculating the shear forces acting on the wick structure. Flow reversal may 
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also affect the entrainment limit. They also found that the pressure profiles 
across the heat pipe were two dimensional. Their study recommended the 
solution of the complete (compressible) two dimensional Navicr-Stokes 
equations to provide an accurate prediction of the pressure drop along the 
pipe. 

McDonald et al. [16] analyzed a vapor flow in an enclosure in the presence 
of air as a noncondensible gas. The coupled governing equations for steady 
state two dimensional flow were numerically solved to study the effects of air 
on the evaporation and condensation rates of the water heat pipe. They con- 
cluded that the two dimensionality of the flow field had a very small effect on 
the concentration of the noncondensible gas and on the temperature distrib- 
utions in the vapor core. 

However, Rohani and Tien [17,18] showed that in liquid metal gas loaded 
heat pipes, a one-dimensional analysis is not adequate. This conclusion was 
explained by averaging the two-dimensional results over the cross section, and 
comparing these with those of one-dimensional analysis [15]. They performed 
a numerical analysis of the steady state two dimensional heat and mass trans- 
fer in the vapor gas region of a gas loaded heat pipe. Their results show that 
the vapor flow in the vicinity of the vapor-gas front is negligibly small and they 
concluded that heat and mass are transfered only by conduction and diffusion 
at the front. But, in an analysis done by Peterson and Tien [19] it was shown 
that natural convection and radial diffusion have a significant effect on the 
noncondensible gas distribution in the vicinity of the vapor-gas front. 
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Gas controlled heat pipes (GCHP) were recently investigated by 
Galaktionov and his colleagues from the Moscow Energy Institute [20.21]. 
In their analysis, they included convection and diffusion across the front under 
heat loads considerably smaller than the sonic limit of heat pipe operation. 
Galaktionov et al. [20] presented a mathematical model for incompressible, 
steady state, two dimensional problem for the heat and mass transfer in the 
region of the vapor-gas front in a GCHP. In their analysis, using boundary 
layer assumptions in the evaporator zone, they applied an integral method and 
approximated the axial velocity by a fourth-order polinomial. 

The same problem was analyzed by Galaktionov and Trukhanova [21], 
taking into account the compressibility of the vapor. Bv defining a stream 
function for the velocity field, the momentum equation was written in terms 
of the vorticity and stream function. The governing equations were solved 
numerically by the Gauss-Seidel method. They also carried out experimental 
studies of a planar gas controlled heat pipe. The results show that the vari- 
ations of temperature in axial and radial directions are quite comparable, 
which clearly demonstrate the two dimensional character of the problem. The 
velocity profiles show the reverse flows in the region of the vapor-gas front. 

Reverse flow in the condenser section was also observed by Ooijen and 
Hoogindoorn [22] who performed experimental and numerical studies on heat 
pipe vapor flow. Their analysis was confined to a steady state, incompressible, 
two dimensional flow. They also found that the total pressure drop over the 
heat pipe, for high radial Reynolds numbers, was higher than that approxi- 
mated by a Poiseuille flow model. 
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Naravana [23] performed a numerical study of incompressible, steady 
state, two dimensional vapor flow. He showed that the pressure terms in the 
Navicr-Stokes equations could be eliminated by integrating the momentum 
equation in the axial direction, on the pipe cross section. In this analysis, using 
boundary layer assumptions, the radial pressure variation was neglected. 

A double-walled concentric heat pipe was proposed by Faghri [24]. Two 
pipes with different diameters create an annular space for the vapor flow. The 
wick structures are attached to the inner surface of the outer pipe and the 
outer surface of the inner pipe. In this kind of heat pipe the area of heat 
transfer into and out of the pipe is increased and its efficiency of performance 
is expected to be higher. Faghri has carried out a numerical study of the vapor 
flow in such a double-walled heat pipe. He used an implicit marching finite 
difference method to solve the incompressible, steady state, two dimensional 
flow problem. His results show that for low Reynolds numbers, viscous effects 
dominate. For radial Reynolds numbers greater than one. pressure decreases 
in the vapor flow direction in the evaporator section, whereas it increases in the 
condenser section due to partial dynamic pressure recovery from the deceler- 
ating flow. For high radial Reynolds numbers, flow reversal flow was shown 
to occur in the condenser section. 

In studies of the dynamic behavior of heat pipes, three different areas of 
concern are startup, shut down and operating transients. Among these the 
startup mode of transient is the most difficult one. As mentioned above, the 
startup of a heat pipe may be encountered with melting of the working fluid 
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at solid state, free molecular flow of vapor and dry out and rewetting of the 
wick. 

Extensive studies of startup and shut down operations of heat pipes 
[25,26,27,28,29] have been carried out at the Los Alamos National Labora- 
tory. The objectives were the applications of high temperature, liquid metal 
heat pipes to space power system heat rejection. Experiments were performed 
on a long heat pipe with a molybdenum container. Using liquid sodium and 
potassium as working fluids, design characteristics of high temperature heat 
pipes have been evaluated. For the temperature range of 700 to 1600 K., the 
useable length of 30 m with length-to-diameter ratio of 300 was found for 
potassium heat pipes, whereas for sodium heal pipes useful lengths of more 
than 40 m with length-to-diameter ratio of 800 were predicted. 

In the numerical studies [29], heat pipe performance was modeled to in- 
clude startup from the frozen state of the working fluid, with provision for free 
molecule flow of the vapor. Vapor core dynamics were simulated by a tran- 
sient one dimensional compressible flow with friction and mass flow from the 
wick. The governing equations were solved using the KACHINA method 
[30]. The numerical code is still under development. In the present status it 
requires 30 minutes of Cray computational time per hour of actual time. 

Bystrov and Goncharov [31] performed experimental and theoretical 
studies on a heat pipe with and without a foreign gas. The steady state oper- 
ation and startup from a frozen state of liquid were analyzed. By averaging 
across the axial cross section, they solved unsteady, one dimensional equations 
for the vapor and liquid phases. Their numerical scheme could be run only for 



very small time steps (I0~ 5 to I0~ 7 sec) due to stability problems and as a re- 
sult needed a large amount of computational time. In order to increase the 
time steps, they used a lumped parameter approach for heat balance and 
avoided calculations of detailed temperature distributions. 

Chang and Colwell [32,33,34] studied the transient operatian of low tem- 
perature heat pipes experimentally and numerically. Experiments were per 
formed with Freon 1 1 as a working fluid to study the dryout and rewetting of 
the liquid phase. In their analysis the vapor core temperature was assumed to 
vary only with time, allowing the solution of a two dimensional conduction 
problem for the liquid and wick structure. 

Ambrose et al. [35] studied pulsed startup of heat pipes. A pulsed startup 
refers to applying a step increase in input heat flux to a heat pipe under steady 
state conditions. They performed experiments on a simple-screen-wicked 
copper-water heat pipe to investigate the transient time and the dryout and 
rewetting behavior. The transient time was found to be dependent on the 
capillary force, the cooling mechanism and the heat capacity of the pipe. 
Dryout of the wick at the evaporator was found to occur only for heat fluxes 
greater than the capillary limit. Dryout and rewetting were modeled using a 
lumped parameter method. Good agreement between the experimental data 
and the theoretical results was reported. 

The works done so far on the transient behavior of heat pipes, mainly con- 
cern startup and shut down operations. However, these modes of transient are 
still not well established. Furthermore, no major work has been done on the 
operating transients heat pipes are expected to undergo. For a comprehensive 
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study of the dynamics of heat pipes, the vapor flow needs to be thoroughly 
investigated both theoretically and experimentally. 

Previous theoretical studies of dynamics of compressible vapor flow have 
considered a transient one-dimensional model. In these studies the friction co- 
efficients of the vapor flow boundaries, have been approximated using a steady 
state two-dimensional model. In addition., steady state two-dimensional ana- 
lyses have shown that one-dimensional models of the vapor flow are not able 
to accurately predict the axial heat and mass transfer and the axial pressure 
drops. Furthermore, two-dimensional steady state studies indicate that flow 
reversal takes place in the condenser section. It is, therefore, important to es- 
tablish the conditions under which this mode occurs during transient oper- 
ations. 
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allows one to induce a liquid-vapor phase transition at low applied heat fluxes. 
A scaling analysis can be used to generalize the results to other fluids. 

Axial and vertical velocity components will be measured using Laser 
Doppler Anemometry. The temperature and concentration profiles in the va- 
por core arc measured simultaneously using twin beam laser holographic 
interferometry. These methods arc described briefly in the following sections. 



Figure 4: Experimental Rectangular Heat Pipe 
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5.2 OPTICAL METHODS 


Optical methods have been used for many years in heat and mass transfer 
measurements due their unique advantages. These methods do not disturb the 
examined process and the information of the whole field is recorded on pho- 
tographs which has advantages over methods of point by point measurements. 
Furthermore, the measurements are inertialess and therefore, very fast phe- 
nomena can be investigated [36]. 

In this section, the use of optical methods in heat pipe experiments arc dis- 
cussed. Laser Doppler anemometry is proposed for measuring the axial and 
vertical velocity components of vapor flow. The advantages of twin beam 
holographic interferometry are discussed as a means of simultaneously meas- 
uring both the temperature and concentration distributions in the vapor core. 


5.2.1 Holographic Interferometry 

Interferometric methods are based on differences in lengths of the optical 
paths. Among these methods, holographic interferometry is the least expensive 
and easiest to construct with a comparative accuracy [36]. However, a highly 
coherent light is needed. The single beam holographic interferometry is first 
illustrated. Then it will be shown how this method can be modified to twin 
beam interferometry. 

Single Beam Interferometry 

Figure 5 is a schematic diagram of a single beam interferometer. The light 
source is a laser. The laser beam is split into an object beam and a reference 
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beam. Each beam is passed through a microscope objective and a collimating 
lens and expanded to a parallel beam. The reference beam falls directly on the 
hologram plate, whereas the object beam passes through the experimental ap- 
paratus and then intersects the reference beam on the hologram plate. 

In real-time holographic interferometry the photographic plate is exposed 
before the experiment starts, and the comparison wave is recorded. The pho- 
tograph plate is developed and fixed while it stays in its place. The comparison 
beam is then reconstructed by illuminating the hologram with the reference 
beam. After the experiment starts a phase change occurs in the object beam 



Figure 5: Schematic of Single Beam Holographic Interferometer 

I. Helium-Neon Laser, 2. Argon Laser, 3. Variable Beam Splitter, 
4. Object Beam, 5. Reference Beam, 6. Microscope, 7. Pinhole, 

8. Collimating Lens. 9. Test Section, 10. Holographic Plate, 

II. Lens, 12. Camera, 13. Air-Suspension Table, 14. Shutter 
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and illumination of this beam on the hologram produces interference fringes. 
Changes in the interference pattern due to changes in the process can be con- 
tinuously observed or photographed. 

Interferometric methods are based on the measurement of variations of the 
refractive index in the test section. These variations may be caused either by 
temperature, concentration or pressure gradient fields. Any of these fields 
alone can then be determined by evaluation of the interferogram. However, 
if the refractive index is changed due to temperature and concentration 
changes simultaneously, an interferogram produced by a single beam 
interferometer may not be adequate to evaluate the temperature and/or con- 
centration fields. This difficulty can be overcome by using two laser beams 
with two different wavelengths. This method is called twin beam (or two- 
wmvclength) holographic interferometry and was first proposed by Mavinger 
and Panknin [37], 

Twin Beam [nter^erometr y 

The first accurate results were obtained by Panknin [38]. Using twin beam 
interferometry he found good results for measurement of the temperature and 
concentration in simultaneous heat and mass transfer along a heated vertical 
plate and a horizontal cylinder with sublimation and the temperature and fuel 
distributions in a flame. 

Figure 6 shows the required experimental arrangement. The twin beam 
interferometer is very similar to that of the single beam shown in Fig. 5. Here, 
two lasers, an Argon laser (Xj = 457.9 nm ) and a Helium-Neon laser 
(X k = 632.8 nm) are used as light sources. The two beams lie on a horizontal 
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plane and arc perpendicular to each other when they intersect on the beam 
splitter. Therefore, both object and reference beams consist of two different 
wavelengths. 

The exposure technique in twin beam interferometry is similar to that of the 
single beam. Before the experiment is started, the comparison wave, which 
now consists of two wavelengths, is recorded on a photograph plate. The 
photograph plate is developed and fixed. Now, if the laser beam with the 
wavelength, say, Xj , is blocked, both reference and object beams will have 
wavelength X k . Then by illuminating the reference beam on the hologram, the 



Figure 6: Schematic of Twin Beam Holographic Interferometer 

I. Helium-Neon Laser, 2. Argon Laser, 3. Variable Beam Splitter, 
4. Object Beam, 5. Reference Beam, 6. Microscope, 7. Pinhole, 

8. Collimating Lens, 9. Test Section, 10. Holographic Plate, 

II. Lens, 12. Camera, 13. Air-Suspension Table, 14. Shutter 
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comparison wave is reconstructed. The object beam which passes through the 
test section is illuminated on the hologram and produces an interference pat- 
tern for wavelength k k . By shining the laser beam of wavelength kj , another 
interference pattern is produced for wavelength kj . 

The main advantages of this technique is that only one plate is developed 
which includes information from both wavelengths, and two interference pat- 
terns are produced by shining each laser beam, separately. Two sets of fringes 
are, then, evaluated to determine the temperature and concentration profiles 
for the tested medium. The evaluation method of fringes is discussed in the 
following section. 


5.2.2 Theory of T win Beam Interferometry 

The following assumptions are made in development of the governing 
equations [38]: 

1. The optical system is perfect, the experimental setup is mechanically 
stable and the lasers are ideal. 

2. The object beams with the wavelengths kj and k k are ideal parallel 
waves. 

3. The variation of the refractive index is only two dimensional and no 
variation in the beam direction. 

4. No reflection of beam due to gradients of the refractive index. 

5. Refractive index is constant outside of the test section and during re- 
cording of the comparison waves. 
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6. The holographic construction is perfect. 

With the above assumptions the interference pattern depends only on the 
variation of the refractive index. In holographic interferometry the object 
beam, passing through the test section at different times, is superimposed on 
the comparison wave, and therefore, reveal the difference in optical 
pathlengths in two exposures. Expressed in multiple 5 of the wavelength X , 
this difference is written as 

S(.x,)vi)-A = /{«(.rjvi)~ n^} (16) 

where l is the width of the test section, in which the refractive index varies be- 
cause of temperature and concentration gradients. The interference pattern 
shows the change of the refractive index between the recorded comparison 
wave (constant temperature and constant concentration C ^ , yield a con- 
stant refractive n x ) and the "measurement" wave. 

An extinction of light occurs for 

\S\ (dark fringes) (17) 

and an amplification for 

|5| = 1,2,3*** (bright fringes) (18) 

The interference fringes are points of the same refractive index change. With 
temperature and concentration gradients the refractive index variation de- 
pends on both gradients. Therefore, a relationship has to be found between 
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the refractive index change due to the temperature gradient and that due to the 
concentration gradient. 

The molar refractivity A 7 (A) , is related to the refractive index by the 
Loren tz- Lorenz equation as follows 


N (A) = 


n (A) 2 — 1 


n(X) 2 + 2 


P 


(19) 


where M is the molecular weight and p is the density. Fortunately the molar 
refractivity does not depend on temperature and pressure, but it varies with 
wavelength. This is the basis for the twin beam interferometry. 

For gases with n ~ 1 , Eq. (19) is approximated quite accurately by the 
Gladstone- Dale equation 

( 20 ) 

For the ideal gases the density, p , can be replaced by the total pressure, p , 
and the temperature, T , to get 

( 21 ) 


where R is the gas constant. The molar refractivity, A r (A) , for a mixture of two 
components is given by 

N{X)=C a N a (X)+C b N b {X) (22) 
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with C a + C b = 1 . In the a bene equation C a and C b are the concentration of 
the components a and b in the mixture and A a (7.) and A T 6 (7.) are the molar 
rcfractivitv of the components in their pure form. 


Combining Eqs. (16), (21), and (22) yields 


S (A- V va) • l = + Crfjr.y) (A'jW - A^))] - -J-A'JA)} (23) 


with C a = 1 - C b , C aoc = 1 and C b<x> = 0 (only component a presents in re- 
cording the comparison wave). For two different wavelengths X } and X k we get 
two equations from which the temperature and concentration distributions will 
be evaluated. These equations are the following 


X i x k 

S (*,}>, A.-) 5 (.xj^;,) = 

^ J AN(Xj) v A N(X k ). 


T(x.y) - 7^ 

3 pi 

N&j) 

N^k) 

T(xy) • T’oo 

2R 

_ A N(Xj) 

AN(X k ) 


X J 

S (X^ylj) 


- S (xy,X k ) 


NM 


C(x.y) 

' 3pl 

AN(Xj) 

AN(X k ) ~ 
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2R 

_ N&j) 
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(25) 


where AA-f/t) = N b {X) - N a {X). 

Figure 7 demonstrates how to evaluate the temperature profile, T{y) , and 
the concentration profile. C(y) , from two interferograms of a boundary layer 
produced by two different wavelengths, X } and X k [38]. For a specified direc- 
tion along the boundary laver, .v = constant , Fig. 7b shows S . S L . S ; • X, and 

w J K J j 
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S k • X k as functions of >• , where 5 y = S(.v = const and 

S k = S( a- = constju k ) . According to Eqs. (24) and (25) the temperature and 
concentration profiles at a location arc found by the following equations 


7lv) - 7^ = 
Qy) ~ ^ 


s r x < 

S t ■ h 

A.Y, 


Sj • Xj 

Sfk 

A\- 

J 

Ajfc 


(26) 


(27) 



a) Interferograms for Xj and X k 

b) Phase lag 5 • X, S 

c) , e) Difference of modified phase lags 

d) , f) Temperature and concentration profiles 

Figure 7: Evaluation of Interferograms [38] 
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where Nj = Njttj) and N k = A ? „(/ k ) . Figures 7c and 7e show the variations of 
5 • a I AN and S • XjN for and A* . Then using these distributions, Eqs. (26) 
and (27) yield the temperature and concentration profiles for that location. 
These profiles are shown in Figs. 7d and 7f. 


5.2.3 Laser Doppler Anemometry 

The laser Doppler anemometry allows the measurement of the local, instanta- 
neous velocity of tracer particles suspended in a flow without disturbing the 
flow. Thus, appropriate particles must exist in the fluid and the relationship 
between the particles and fluid velocity must be known. The task of a laser 
Doppler anemometer falls into three parts. 

• Generating two coherent laser beams which interfere in a probe volume 

• Detecting the light intensity scattered by tracer particles crossing the 
measuring control volume 

• Analysing the Doppler frequency and calculating the velocity of the 
particles 

Figure 8 shows the required components of a laser Doppler anemometer and 
the following is a brief description of the basic principles of laser Doppler 
anemometry which are carefully outlined in Durst et al. [39]. 

Two coherent laser beams with wavelength X and a diameter d form an 
ellipsoidal volume at the beam intersection point (see Figs. 9 and 10). The size 
of the measuring control volume is given by 
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5 X = d ; 


d 

sin (f> 


(28) 




cos <fi 


where <p is the angle between the two beams and d x , d y and <3. arc the axes of 
the ellipsoid. Figure 10 shows the principles of the "interference fringe" model 
which was proposed by Rudd [40-]. 

Two coherent laser beams having plane wave fronts intersect at an angle, 
<£. This yields a pattern of plane interference fringes as shown in Figure 10. 
The fringe spacing , A z, is proportional to the wavelength, X , of the light and 
inversely proportional to the half angle between two beams 



Figure 8: LAD Optics from DANTEC Manual 

1. Beam Splitter 2. Bragg Cell Unit 3. Beam Displacer 
4. Beam Expander 5. Front Lens 6. Pinhole Unit 
7. Photomultiplier Optics 8. Photomultiplier 
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As a particle moves through the region of interference fringes it will scatter 
light whose intensity will vary according to the light intensity variation inside 
the bisector of the two beams. These light intensity variations have the 
Doppler frequency given by 

U 2Usin (<f>l2) 


where U is the particle velocity component perpendicular to the fringes. 

The scattered light can be detected by a photomultiplier and the velocity 
component normal to the interference fringes can be determined by analyzing 
the Doppler frequency. This frequency is ambiguous because it docs not indi- 
cate in which direction a particle moves. This lack of information can yield 
large errors if the flow direction is not known (e.g. in ease of recirculating flow) 



Figure 9: Measuring Control Volume 



or if the mean velocity is near-zero and the fluctuations arc larger than the 
mean value in turbulent flow. 

Frequency shifting allows the 180 degree ambiguity in the direction of a 
measured velocity to be resolved. The frequency shift is obtained when the two 
intersecting laser-beams have different frequencies. The pre>,cnce of a fre- 
quency difference, v 5 , results in a movement of the fringe pattern with velocity 

U* 

U$= — v^Az (31) 

A frequency shift, applied to one of the laser beams, will result in an increase 
or decrease in the measured velocity of a particle depending on the direction 
of the movement. In this ease the velocity can be determined clearly from the 
detected Doppler frequency, . If the velocity of the fringes, caused by the 
shift frequency, is greater than the particle velocity 


Z 

t 



Figure 10: Interference Fringes in the Control Volume 
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/ 


(32) 


t' = (V0-V S )- 


2 sin {cp! 2) 


A lower Dopplcr-frequency than the shift frequency indicates that a particle 
moves in the same direction as the fringe pattern. If a particle moves in the 
opposite direction, the Doppler frequency will be larger than the shift fre- 
quency. 

Laser Doppler anemometry depends on signals from particles suspended in 
the flow. Thus, it is important to know how well the motion of the particles 
represent the motion of the fluid. In our case there are two reasons for velocity 
differences: 

• The particle terminal velocity is different from the vapor velocity be- 
cause of greater gravitational forces on the particle as compared to the 
vapor. 

• Drag forces caused by accelerations and velocity fluctuations can result 
in particle/vapor mismatch. 

According to Tchen [41], the terminal velocity, C f , of a particle in a fluid with 
the constant velocity U f parallel to the gravity g is 

2 

Up — Uf— Cj with C, = f-£-(p p -pj) (33) 


The velocity difference C, between the particle U p and the fluid Uf depends on 
the density difference, p p — p f , the viscosity of the fluid, pp and the squared 
radius of the particles rj . For oil droplets (r p ^\pm) in gas, at standard state. 
C.~ 1 .5 • 10 -4 m s. 
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According to Dring [42] the velocity difference between a particle and fluid 
caused by high-frequency velocity oscillations, depends on the Stokes number 


P p u4r p 
18 ft/ 


(34) 


For Stokes number St < 0.14 the particles follow an oscillation with the fre- 
quency co with less than 1% velocity difference. For oil droplets in a gas. at 
standard states, 12 kHz. 


5.3 PRELIMINARY RESULTS 

Experiments have been carried out to visualize the vapor flow patterns in 
evaporation, adiabatic and condensation regions of the heat pipe after a sud- 
den increase in input heat flux has been applied. A boundary layer of vapor 
is established very quickly on top of the wick surface all over the heat pipe. 
The thickness of this layer decreases from the evaporation section to the 
condensation region. The vapor boundary layer grows quickly and fills the 
whole space. Depending on the input heat flux, the whole process takes place 
in a few seconds. For low input heat flux, the vapor flow is found to be smooth 
and laminar. For higher heat fluxes the vapor boundary layer grows faster. If 
the input heat was high enough a few regular vortices were observed on top 
of the evaporation part and they expanded to the entrance region of the 
adiabatic part. The vapor flow in the adiabatic and condensation regions un- 
der a high heat flux arc still laminar but reverse flow on top of the condenser 
is suspected. 
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Transient processes in the vapor core of heat pipes are very fast (few sec- 
onds). Therefore, the mentioned optical methods which provide inertialess 
measurements, are the proper tools to obtain accurate results. The temper- 
ature and concentration profiles are evaluated from the holograms developed 
by the twin beam interferometry. In this method, the temperature is evaluated 
from the difference between the phase shifts corresponding to two different 
wavelengths (see Eq. (26)). This difference is usually very small, therefore, the 
two wavelengths used should be as far apart as possible. 

In laser Doppler anemometry, seeding of particles with proper density and 
size is a challenging problem. In systems like heat pipes, in which the working 
fluid undergoes a phase change, it is not possible to have particles recirculating 
with the fluid flow. In this case, particles can be injected carefully into the 
vapor flow in the evaporation end and collected from the condensation end. 
But the disadvantage of this technique is that the measured velocities at both 
ends are not accurate. 

When the evaporation and condensation rates are high, small liquid drop- 
lets are present in the vapor flow. These droplets flow with the vapor and are 
able to scatter the laser beam. Preliminary velocity profiles in both the 
evaporator and condenser have been obtained to demonstrate that the method 
will be effective. Therefore, no seeding will be required in this case. 
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Chapter VI 

THEORETICAL MODEL OF THE VAPOR CORE 


Flow patterns found in a heat pipe are shown schematically in Fig. 11. The 
working fluid evaporates in the evaporation zone. The local pressure increases 
and the vapor flows towards the condensation zone, at which it condenses back 
to the liquid phase. The liquid flows from the condenser to the evaporator by 
capillary force through the wick structure. The main concern of this analysis 
is the vapor core response to changes in the evaporation and condensation 
rates due to a sudden increase, or decrease, in the input heat flux, or the 
condenser temperature. 

In this section the dynamic behavior of the vapor flow is analysed using a 
transient two dimensional model. The proper equations which govern the 
model problem will be written. A numerical scheme will be proposed to solve 
the governing equations. Then, presenting some preliminary results of the 
numerical calculations, the continuation of the analysis to more rather com- 
plete solution of the problem, will be proposed. 
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Figure 1 1 : Flow Patterns in a Heat Pipe 
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6.1 GOVERNING EQUATIONS 

The vapor How in the heat pipe vapor core, shown in Fig. 1 1, is modeled as a 
channel How shown in Fig. 12. The circular boundary of the channel is a thin 
porous medium which contains the liquid. The input heat flux to the 
evaporator and the temperature of the outer surface of the condenser are 
specified. The planar side wails are assumed adiabatic. 

The equations governing the vapor flow arc time dependent, viscous, 
compressible momentum, continuity and energy equations. An equation of 
state is used to relate pressure, density and temperature within the vapor core. 
These equations in axisymmetric coordinates (r,x) are 



Figure 12: Vapor Core Model of a Heat Pipe 
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6 . 1.1 Boundary Conditions 

The boundaries of ihc vapor core are shown in Fig. 12. The nonslip condition 
for velocity and adiabatic condition for temperature are assumed on the side 
walls, i.e., 

@ .t = 0 and x = L . 

u = 0 , V = 0 , -P- = 0 (40) 

ox 


On the center line the symmetry condition implies. 


@ r = 0 , 


du 

dr 


= 0 , v = 0 , 


dr 

dr 


= 0 


(41) 


The boundary conditions on the liquid-vapor interface are the challenging 
ones. The liquid flow is assumed in a porous medium with thickness S, which 
is much smaller than the vapor core radius r 0 . The axial velocity is assumed 
zero on this boundary. That is, 


'S r = 'o 


u = 0 


(42) 


In order to assign boundary conditions for the temperature and vertical veloc- 
ity, the liquid-vapor interface is divided into three regions. In the evaporation 
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zone the input heat flux, q, is a given parameter and the input flow is approx- 
imated as, 

@ r — r 0 and 0 < x < L e , 


pv = m ~ 


<? 


(43) 


In the above equation, h fg is the heat of evaporization and conduction in the 
liquid layer is ignored. The temperature is assumed to be the saturation tem- 
perature of the liquid corresponding to the interface pressure. That is, 

@ r — r 0 and 0 < jc < L e , 

T = T sakP) ( 44 ) 


In the adiabatic zone the boundary conditions are 


@ r = r 0 and L e < x < (L e + L c ), 

v = 0 , -^=0 (45) 

dr 

In the condensation zone the temperature of the outer surface, T c , is a given 
parameter. By equating the heat of evaporation to the heat conduction in the 
liquid layer, the outflow from this region is approximated as, 

r = r {) and (L e + L c ) < x < L. 
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px = m — 


(46) 


k eff( T-T c \ 

b \ h f^ J 


T= T sa £p) (47) 

where, k eff is the effective conductivity of the liquid layer and the temperature 
is assumed to be the saturation temperature corresponding to the pressure at 
the interface. 


6.2 MIXED GALERKIN-FINITE DIFFERENCE METHOD 

A mixed Galerkin-finite difference technique was used for the solution of the 
steady state equations, see Ref. [44]. In each partial differential equation the 
equation of state was substituted for the pressure, eliminating the pressure as 
a variable. This gave a system of four equations in four unknowns. Galerkin 
method was then applied, resulting in a pseudo two dimensional set of 
equations. Application of Galerkin method [45] involved choosing trial func- 
tions to represent each of the variables and substituting these in the equations. 
For the correct implementation of the Galerkin method, the trial functions 
chosen must have three requirements met. 

1. They must be linearly independent. 

2. They are part of a complete set of functions. 

3. They must satisfy the boundary conditions exactly. 
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In this case the trial functions described the radial variation of the variables 
and exactly matched the boundary conditions in the radial direction. The trial 
functions used were as follows: 

p(x,r,x)=p w + 'Y}A x ' x M r ) 

n 

\ix,rj)=f- i (r)V w + y^v^i^r) 
n 

u{x,rx) = '^u r lx,xy l (r) 
n 

T{x,r,T)= T w +y^T r lx,xy ] (r) 
n 

where the basis functions f l (r),f 2 (r) and f 3 (r) are defined as: 

f\( r ) = cos^(2n — 1)-^-) 

/ 3 (r)= (ir ) 2 

Note that the functions f n are chosen to satisfy the radial boundary conditions. 
Three of the trial functions include a term in front of the summation of coef- 
ficients. This additional term, denoted with the subscript u\ represents the 
known value of the variable at the wall, where r= /?,. At the wall the radial 
functionality of the trial function becomes zero, so that the variable is com- 
pletely defined by the wall term. A function of radius,^ was included with the 
wall term for radial velocity. This was to insure that the radial velocity went 
to zero at the centerline. 
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Substitution of the trial functions into the differential equations creates the 
residual. R. Finally, weighting functions. F w . are chosen for each equation and 
the inner product of each residual and its weighting function is calculated. 
Weighting functions chosen for this analysis are given below: 

Mass 

F wl =fi(r)=cos((2k - 1 )-^-) 

X-Momentum 

F w 2 = f\( r ) = cos((2/r - 1)^-) 

R-Momentum 

F w3 =M r ) = sin^kn 

Energy 

F w4 = f\( r ) = cos{(2k - 1 ) 

Ideally the weighting functions should be orthogonal with the trial functions, 
but in the problem at hand that was not possible. Here, we define the inner 
product as: 

(R,F W )= f F w Rrdr 
•V 

where R is the residual, F w is the weighting function and the product is inte- 
grated over the cross section. These calculations arc carried out in detail in 
Reference 44. A set of nonlinear differential equations (with respect to the 
axial direction) result that can be solved using any of several methods. 
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6.2.1 Results 


Although more promising than the finite difference method described in sec- 
tion 6.3, a number of difficulties were encountered and the press of time con- 
straints did not allow them to be overcome. It is still our view that the mixed 
finite difference-Galerkin method or full Galerkin method is the only method 
that will yield the needed computational speed. The following paragraphs 
document some of our struggles to "make it work". 

Converged results from the analysis have not provided any more insight 
into heat pipe vapor dynamics than was previously known. Several types of 
vapor flow in the heat pipe have been input to the code to see if accurate and 
meaningful results could be obtained. Solutions used for checkout corre- 
sponded to a 'zero solution', a stagnant tube, fully developed pipe flow and 
finally the heat pipe operating mode. 

Corresponding to a zero solution' are Dirichlet boundary conditions of zero 
for density, radial velocity and temperature at the wall. The corresponding 
internal coefficients for the Galerkin representation should then be calculated 
to be zero as well. A stagnant tube corresponds to the geometry of a sealed 
tube at constant, uniform temperature and density. Hence, values for density 
and temperature through the pipe are equal to the wall value. Coefficients for 
the Galerkin representation of these variables should be zero. Consistent with 
the stagnant tube, all vapor velocities should also be zero. 

Determining boundary conditions consistent with flow in a pipe required 
that the program be modified slightly. Modifications were primarily in regards 
to the axial boundary conditions on axial velocity. Previous conditions stated 
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that all velocity components at x = 0 and x = L were zero due to the solid 
boundary at the ends. For pipe flow, it was intended to model a small section 
cut from a long, straight pipe. Thus, in the model, a fully developed for the 
axial velocity could be used at x = 0, and at first was also used at x = L. 
However, the conditions proved to be inconsistent because vapor flow was be- 
ing forced through a tube, but with constant density and temperature there 
was no pressure drop. A remedy to this problem was to use a Neumann 

boundary condition at the right end of the pipe, = 0 at x = L , and set ei- 

dx 

thcr density or temperature to be decreasing axially from x = 0 to x = L, with 
a linear profile. In this manner, the actual value of the axial velocity was not 
set, so it could adjust to the changing pressure along the pipe. 

Finally the heat pipe operation was implemented. Under this mode of op- 
eration, the axial velocity is set equal to zero at all the boundaries while the 
radial velocity is set equal to zero at the two ends, x = 0 and x = L. In the 
evaporator the radial velocity is set equal to a finite negative value at the wall, 
and to a finite positive value at the wall in the condenser. At the wall in the 
adiabatic section a zero radial velocity was assumed corresponding to no mass 
evaporation or condensation. Because this was a steady state analysis, the 
total mass flow in the evaporator was equal to the total mass flow in the 
condenser. Density at the wall was held constant for this try and the wall 
temperature was given a profile which decreased axially toward the condenser. 
This was consistent with a heat pipe operating under a steady heat flux and 
allowed for decreasing pressure in the direction of flow. 
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An instance may be documented in which the current equation set will fail 
as it is. For the zero solution', i.e.. all boundary conditions are zero, so the 
coefficients of the Galcrkin trial functions are also zero, the equation vector 
on the right hand side of the iteration equation is indeed zero, but for a sol- 
ution, the Jacobian must be calculated and non-zero terms must be obtained. 
Jacobian elements for all terms on the diagonal may be determined for the 
energy equation and the two momentum equations. However, the conserva- 
tion of mass equation will have no non-zero terms, because all Jacobian de- 
rivatives of this equation still have a coefficient of wall term which is zero. 

The code was run for the 'zero solution'. Two sets of initial conditions were 
used for the Galerkin coefficients. As a first guess the coefficients were initially 
set to zero. This proved to exactly solve the set of equations, as only one iter- 
ation w’as used and the zero solution was returned. Then the initial guess for 
the coefficients was set to 0.0005. With this guess the program required two 
iterations to converge with a maximum iteration error of 6.5 • 10 -5 on the 
temperature. In this case, however, the coefficients were non-zero, with the 
largest magnitude being 5.0 • 10 -4 for the density, 3.1 • 10 -8 for the axial ve- 
locity, 1.0 • 10 -8 for the radial velocity and 8.4 • 10 -4 for the temperature. 

Solution for a stagnant tube required 21 iterations from an initial guess of 
0.0005. The maximum iteration error was 9.6 • 10~ 5 . Again, this obtained for 
the temperature. The maximum magnitude of the converged values was 
1.5 • 10~ J for the density. 2.5 • 10 -7 for the axial velocity, 1.7 • 10 -6 for the ra- 
dial velocity and 4.7 • 10~ 4 for the temperature. For comparison, wall values 
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for density and temperature were 1.0 and 3.0. respectively, and the calculated, 
normalized sonic velocity was 2.8. 

These two cases are the onlv ones the code would definitely converse on. 
They were obtained on an eleven grid matrix with a spacing of 0.01. Following 
these two cases, the pipe flow configuration was attempted. Using the same 
grid network with Dirichlet boundary conditions at x = 0 and Neumann con- 
ditions at x = L for the axial velocity, no solution was obtained. In this con- 
figuration only the axial velocity coefficients were updated. Radial velocity was 
everywhere set equal to zero as were the coefficients of the density and tem- 
perature. Hence, the density and temperature were axially and radially con- 
stant. Iterations on the axial velocity diverged slowly at first, then rapidly 
around the 18th iteration. 

An additional pipe flow case was considered in which both axial and radial 
velocities were updated at each iteration. Using the same grid network the 
solution diverged again, this time after 32 iterations. 

Finally, the code was run for the heat pipe configuration, using constant 
density and injected vapor velocity of -0.01. Condensing vapor in the 
condenser was given a velocity of +0.01. Wall temperature varied axially 
from 3.2 in the evaporator to 2.S in the condenser. Updating all variables and 
using an initial guess of of 0.0 for all coefficients, the solutions diverged after 
9 iterations. A damping factor of 0.77 was included. Here, the calculated co- 
efficients after the 9th iteration displayed a noticeable amount of cell Reynolds 
number problems. To try and correct the problem, the grid spacing was cut in 
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half. Using the same boundary conditions and the initial guess, the solution 
appeared to be oscillating after 40 iterations and failed to converge. 

In general, results from using the code for the heat pipe configuration have 
all failed. One problem evident from the output has been a cell Reynolds 
problem. However, increasing the number of nodes, although permitting more 
iterations without diverging, still does not allow consistent, meaningful results 
to be calculated. 

Although grid spacing and the resulting cell Reynolds problems may have 
been a problem, it does not appear to be the only one. Some equation formu- 
lations have instabilities within them and one way to arrive at consistent sol- 
utions is by including damping. A method of upwind differencing was utilized. 
This too seemed to make little difference to the computed solutions. 

Completeness of the basis functions used to represent the variables in the 
Galerkin expansion is another point in question. For the cylindrical coordinate 
system it probably would have been better to use Bessel functions as the basis 
functions. 

It looks like the boundary conditions themselves may have caused problem 
which the mixed Galerkin-finite difference method does not handle well. 
Boundary conditions on the radial velocity are sharply discontinuous, with the 
velocity being negative in the evaporator, zero in the adiabatic region and 
positive in the condenser. Changes between the regions occur as sharp jumps. 
Trial functions chosen to represent this variable could not satisfy these condi- 
tions with a single basis function alone, but relied on the constant term added 
to this basis function to match the requirements completely. However, con- 
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stant term still poses sharp discontinuities in the axial direction. It may be that 
the failure of the basis functions alone to satisfy the boundary conditions 
causes the Galerkin method to not converge to a solution. It appears the 
Galerkin method is not well suited to handle the boundary conditions associ- 
ated with this type of a compressible flow problem. 

In conclusion, there are several aspects of this work that could have been 
valuable, had a solution be obtained. The basic problem studied, compressible 
heat pipe vapor dynamics in two dimensions, is a topic of current interest in 
many communities around the world and a capability of detailed analysis is 
needed. Secondly, the method of solution, mixed Galerkin-finite difference 
method, is not one that is often applied to a compressible flow problem. Suc- 
cessful application of this method as well as the direct solution of the steady 
state equations could have been a contribution of interest to the field of nu- 
merical fluid mechanics. 


6.3 PROPOSED SOLUTION METHOD 

Equations (35) to (39) are five equations in five unknowns, namely, p , u , v , 
T and p . The boundary conditions associated with these equations are Eqs. 
(40) to (47). The boundary conditions on the vapor-liquid interface depend on 
the interface pressure and pressure is not available on the boundaries. How- 
ever, in the numerical method described below, the pressure does not need to 
be specified on the boundaries. Then, after completion of the calculations, an 
interpolation method is used to determine the pressure on the boundaries. 




X 


— locations for u, f locations for v, 
o locations for p, T and p 

Figure 13: Staggered Grid 

The numerical scheme used here is based on the SIMPLE method described 
by Patankar [43]. In this method all the variables are not evaluated on the 
same point, but a staggered grid, as shown in Fig. 13, is used. By choosing 
control volumes which have the cell centers at their center, the velocity com- 
ponents are evaluated for the points that lie on the faces of the control vol- 
umes. The locations for u and v are shown by short arrows in the Figure. The 
other dependent variables, namely, p , T and p , are calculated for the cell 
centers shown by small circles. The boundaries do not lie on the cell centers. 
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The algorithm for calculations is the following : 

The system pressure is assumed to be composed of two parts 

p{x,r,t) = p{t) + p(x,r y t) (48) 

where p(t) is the space averaged system pressure whereas p(x,r.t) includes the 
space variation of the pressure. Assuming a pressure distribution, the mo- 
mentum equations (Eqs. (36), (37)) are solve for velocity components u and v 
at a time step. Then using the continuity equation and density p from the last 
iteration on temperature, the pressure is corrected to be used again in mo- 
mentum equations. Upon accomplishing these iterations, the velocity field and 
p{x,r,t) will be found. 

Using the velocity field, the energy equation (Eq. (38)) is then solved for 
temperature. In order to calculate density from the equation of state, and to 
evaluate the boundary conditions, we need the absolute value of the pressure. 
Global mass balance is used to calculate the absolute pressure. The total mass, 
for an ideal gas, m , is related to the pressure by 



where V is the volume. The increase in the total mass in one time step is 


m{t 4- At) - m(t) = (m jn - rh ou[ )At 


(50) 
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where m in and m ou . arc the input and output mass fluxes, respectively, and arc 
calculated from the boundary conditions. Equations (49) and (50) are com- 
bined to yield 


{m in - m ou .)At = 


{— 

\RT I 


t+At 


t+At 


-Eli-V 

RT\,J 


(51) 


Thus Eq. (51) is used to evaluate the absolute pressure of the system. For a 
real gas, using an equation of state, the above procedure is followed to calcu- 
late the total pressure. 


6.4 PRELIMINARY RESULTS 

The computational procedure described above was implemented for a two di- 
mensional time dependent compressible flow with the same type of heat pipe 
boundary conditions as those mentioned above. The computations were per- 
formed for several evaporator heat fluxes with water as the working fluid and 
pipe dimensions of r 0 = 2.5 cm and L = 3 r 0 . The computational grid spacings 
were A x/r 0 = 0.1 and Ay I r 0 = 0.05. The temperature and pressure variations, 
and consequently, the density variation, were found to be quite small inside the 
pipe. This is expected of systems in which condensation takes place without a 
non-condensible gas. 

The steady state results are shown in Figs. 14. 15 and 16. The steady state 
was assumed to be reached when the relative change of all variables (error 
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Figure 1 5: Schematic of Flow Streamlines 

(a) low input heat flux, (b) high input heat flux. 

norm) within a time step was less than 0.5%. Figures 14 show the axial veloc- 
ity profiles at different axial cross sections along the pipe, for low input heat 
flux, Fig. 14(a), and for high input heat flux, Fig. 14(b). The flow stream lines 
for these two cases arc shown schematically in Fig. 15. 


60 


y/b=0 



Figure 16: Vertical Velocity Profiles 

When the input heat flux is relatively low, the evaporation rate is low and 
the vapor flows smoothly towards the condensation region (Figs. 14(a) and 
15(a)). However, with high input heat flux the evaporation rate is high, the 
vapor is ejected with a high momentum from the evaporator and under steady 
state conditions, condenses at a high rate at the condenser. The high momen- 
tum flow at two ends of the adiabatic section causes a circulation flow in this 
region. That is, on the adiabatic surface there is a reverse flow. In addition, 
in the condensation region the high momentum vapor flow impacts the rigid 
end wall, returns and then condenses in the liquid region. This return flow is 
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0- initial profile (t = 0), 

I - r = 10- 4 , 

2- i = 2.5* 10- 3 , 

3- i = 10" 2 , 

4- steady state (r = 2.5* 10 -2 ) 
where r = tajb 1 . 

Figure 17: Axial Velocity Profile at Different Time Steps 

shown in Figs. 14(b) and 15(b). The vertical velocity profiles for different ver- 


tical cross sections are shown in Fig 16. The flow is clearly two dimensional. 


Transient behavior of vapor flow was studied by using a steady state sol- 


ution as the initial conditions for different input parameters. As an example. 


Fig. 17 shows the axial velocity profiles at the middle of the pipe, xjL = 0.5 , 


at different time steps. Curve 0 on this Figure is the steady state solution for 


Re= 100 , where Reynolds number is defined as 
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m in r 0 

Re= M 

and m in is the input mass flux. Curves 1 to 4 show this velocity profile at dif- 
ferent time steps after a sudden increase in input heat flux for Re = 200. At 
the outset, the increase in input heat flux overcomes the reverse flow, see curve 

1 . As the velocity profile develops with time, the steady state profile, curve 4, 
shows a higher reverse flow at the liquid-vapor interface. 


6.5 PROPOSED ANALYSIS 

The dynamic behavior of vapor flow in heat pipes is studied for transient op- 
erations. The above analysis needs to be modified as follows: 

1. The existing computer code is able to handle two dimensional transient 
vapor flow in Cartesian coordinates. The code must be modified to 
axisymmetric coordinates for analyzing circular heat pipes. 

2. Interesting phenomena have been detected in the adiabatic and 
condensation regions. In these regions high input heat flux may cause 
a reverse flow. The reverse flow results in negative shear force on the 
wick structure. A parametric study is needed to explore the reverse flow 
versus design parameters like the evaporator, condenser, adiabatic and 
total length of the heat pipe for different input heat fluxes and 
condenser temperatures. 

3. Effects of the presence of a non-condensible gas in the vapor core on the 
condensation rate and the heat pipe operation will be also investigated. 
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4. Many heat pipe applications result in irregular geometries. Grid gen- 
eration techniques need to be incorporated into the code to accommo- 
date them. 

5. The vapor flow has been dealt with using simplified boundary condi- 
tions. The liquid phase model, see Chapter VII, needs to be used to 
generate the proper boundary conditions. 
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Chapter VII 

THEORETICAL MODEL OF THE LIQUID PHASE 

As previously discussed, most research in heat pipe dynamics has focused 
on the startup transient. An experimental investigation has also been per- 
formed for the shutdown transient. Only modest efforts have been expended 
to date that consider operational transients. These efforts generally address 
two issues: the capability to achieve nominal operation under specified condi- 
tions, and the resulting axial temperature distribution. The processes are as- 
sumed to occur quasisteadily such that time dependent characteristics are not 
considered. Despite these research efforts, a comprehensive analytical model 
does not exist that captures the physical processes while being economical 
enough for required computations. 

The objective of this part of the study is to develop an analytical model of 
the liquid phase of a high temperature heat pipe. The model is intended to be 
coupled to a vapor phase model for the complete solution of the heat pipe 
problem. The mathematical equations arc formulated consistent with physical 
processes while allowing a computationally efficient solution. The model sim- 
ulates time dependent characteristics of concern to the liquid phase, including: 
input, phase change, and output heat fluxes, liquid temperatures, container 
temperatures, liquid velocities, and liquid pressures. 
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This research considers requirements for high temperature, high thermal 
flux operation with heat transport over long distances and small temperature 
drops. These requirements indicate the need for a heat pipe with a high per- 
formance composite wick. A composite wick consists of small pores at the 
liquid-vapor interface and large pores in the direction of bulk liquid flow. 
Recall that the maximum capillary pumping capability is limited by the largest 
capillary surface at the liquid-vapor interface as shown by Eq. (4). The small 
pores at the liquid-vapor interface establish a high capillary pumping capabil- 
ity that drives the circulation of working fluid. The larger pores in the direc- 
tion of flow provide a low resistance flow path for the liquid. 

The heat pipe selected for this study uses an annular wick configuration as 
shown in Fig. 18. The annular wick is an excellent configuration for liquid 
metals, is popular for experimental work, and involves geometry easily mod- 
eled in cylindrical coordinates [27,28]. The wick configuration consists of se- 
veral layers of fine pore screen pressed together and concentrically installed in 
the pipe to form an annular flow channel as shown in Fig. 19. The annulus 
forms the low resistance flow channel while the fine pore screen tube forms the 
high pumping capability boundary between the liquid and vapor regions. 

The capillary wick introduces some analytical although not conceptual dif- 
ficulty. A thin screen tube consisting of several layers of pressed screen forms 
a tortuous flow channel with no slip conditions on the surface of the screen 
wires. Explicit analysis of the flow including consideration of the large number 
of wire surfaces involved would be difficult if at all practical considering the 
goal of performing transient analysis. A macroscopic approach such as 
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Figure 18: Heat Pipe with High Performance Composite Annular Wick 



Figure 19: Cross Section of High Performance Annular Wick Configuration 
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Darcy's Law for flow in a porous medium could be used to alleviate the com- 
plexity. However the axial flow in the pressed screen is expected to be small 
relative to axial flow in the open annular channel. Indeed, the very purpose 
of the annular composite wick configuration is to provide a low resistance flow 
channel in the annulus for circulation of liquid. The screen tube is used only 
to provide a fine capillary surface at the liquid-vapor interface. The screen 
tube can then be approximated as a single layer of screen with a characteristic 
pore size. Flow normal to the screen is restricted by a characteristic specified 
porosity, while tangential flow is governed by alternating free slip and no slip 
conditions. The free slip condition provides the tangential shear communi- 
cation between the liquid and vapor phases. The wick is approximated then 
as an infinitesimally thin porous surface. 


7.1 THE LIQUID PHASE MATHEMATICAL FORMULATION 

7.1.1 A ssumptions 

The governing equations to be presented below will be reduced by recog- 
nizing unique characteristics of the system to be analyzed. The system is as- 
sumed to operate in an essentially gravity free environment. Other body forces 
such as those due to translational and rotational acceleration will also be neg- 
lected. 

The heat pipe is assumed to be installed with azimuthal symmetry in the 
particular application. The thermal energy source adds heat uniformly around 
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the evaporator circumference while the heat sink extracts heat uniformly 
around the condenser circumference. The resulting uniform conditions in the 
azimuthal direction imply there are no velocity components or gradients in the 
azimuthal direction. 

Liquid flow in a heat pipe wick is generally low speed, laminar, and 
incompressible. Since the liquid flow field is relatively simple and the cross- 
stream dimension of the liquid region is much smaller than the strearmvise di- 
mension, the cross-stream (i.e., radial) pressure gradient is assumed to be 
negligible. Also since the flow is low speed and incompressible, viscous energy 
dissipation can be neglected. 

In summary, the assumptions to be used for reducing the governing 
equations are: 

1. Negligible body forces: f r =f x = Q 

2. Azimuthal symmetry 

a. No azimuthal flow: w = 0 

3 * 

b. No azimuthal gradients: — — = 0 

36 

3P 

3. Negligible radial pressure gradient: — — = 0 

dr 

4. Negligible viscous dissipation 

dp, 

5. Incompressible liquid: — — • = 0 

6. Laminar liquid flow 

Thermoplmical properties arc considered to be uniform functions evaluated 
at a reference temperature. The reference temperature is assumed to change 
slowly with time such that properties can be treated as constants in the gov- 


erning equations. 
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7.1.2 Governing Equations 

Temperatures of the heat pipe container and the temperatures, pressures, 
and velocities of the liquid phase are required to describe the dynamic behavior 
of the liquid phase of a high temperature heat pipe under transient operating 
conditions. The governing equations describing the behavior of these param- 
eters consist of the conservation of mass, momentum, and energy. 

The full condition governing conservation of mass of an incompressible fluid 
is 


_| iL + 4 L + J.i^ + v =fl (52 ) 

ox or r o<p r 


where u is the velocity in the axial (*) direction, v is the radial (r) velocity, and 
w is the azimuthal direction ($) velocity. Imposing the assumption of 
azimuthal symmetry reduces Eq. (52) to the familiar two dimensional form 


5u 

dx 



(53) 


Full conservation of radial momentum for an incompressible fluid with con- 
stant properties is given by 


- 2 , - 

ov duv gv J_ gwv 

dt dx dr r d4> 



1 dP 

P 5r 


+ v 



+ 


g~v 1_ 

-2 2 = ,2 
Or r d(p 


1 dv 
r dr 


2 dw 

r 2 c4> 



(54) 
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where t is the time, p is density. P is pressure, and v is kinematic viscosity. 
By assuming negligible radial pressure gradient, the conservation of radial 
momentum equation is replaced by 



(55) 


Under this approximation the radial momentum equation is not used to define 
liquid dynamics so that pressure is determined from the continuity equation 
and the axial momentum equation. 

Conservation of axial momentum for an incompressible fluid with constant 
properties is given by 


du cu~_ dux 1 duw uv _ 

dt dx dr r d<p r 

__ 1 dP f d\t d~u 1 dju 1 cu\ 

P dx \<3x 2 dr 2 r 2 d(f> 2 r dr ) 


+ fx 


(56) 


Assuming negligible body forces, azimuthal symmetry and no radial pressure 
gradient reduces Eq. (56) to 


du d\P_ duv uv_ 

dt dx dr r 


1 dP 

P dx 



+ 




(57) 


Full conservation of energy for an incompressible, constant property fluid 
is given by 
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Neglecting viscous dissipation and assuming azimuthal symmetry 


cT &T dT , 

— h u— h v-^— = a 

ct cx or 
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or 


(59) 


The conservation of mass and momentum equations are solved subject to 
the boundary conditions using an integral method, whereas the conservation 
of energy equation is solved using the finite difference method. The computa- 
tional implementation of the solution approach will be discussed in section 7.3. 


7.1.3 Boundary Conditions 

The heat pipe boundary conditions can be considered as either external or 
internal. The external conditions refer to the communication of the heat pipe 
external surface with the surrounding environment. The internal conditions 
are of two types: fluid-solid and fluid-fluid. The fluid-solid conditions refer to 
the interaction between the liquid and either the pipe internal surface or the 
wick. Fluid-fluid conditions govern the interaction between the liquid and 
vapor phases. 

The heat pipe external boundary can be divided into three convenient re- 
gions as shown in Fig. 18. The evaporator region absorbs heat from the envi- 



ronmcnt, the adiabatic region exchanges no heat with the environment, and 
the condenser rejects heat to the environment. The evaporator surface is as- 
sumed to be exposed to a known surface heat flux condition, while the 
condenser exchanges heat through radiation 1 to a sink of known temperature 
[26]. The conditions on the pipe external wall are then 


AXIAL RANGE RADIAL RANGE 

0 < A- < L e r = R w 


CONDITION 


-*4r- = Rinix, 0 

Or 


L e < c x < L e + L a r 


8T 

dr 


= 0 


L e + L a < x < L r = R w 



Rrad 


where q in is specified and q rad is to be determined from q rad = — T \ ) as- 
suming environmental sink temperature T e is known. 


1 The condenser could generally exchange heat with the environment 
through conduction, radiation, and convection. However, the study of 
radiation loaded elements is of current research interest for space based 
applications [25]. 
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Long heat pipes with large evaporator and condenser circumferential sur- 
face areas relative to the surface area of pipe endcaps will normally have neg- 
ligible effects due to conduction in the end caps. With this assumption, the 
boundary conditions are: 


AXIAL RANGE 

RADIAL RANGE 

CONDITION 

* 

II 

o 

Ri<r< R w 

* H 

ii 

o 

x = L 

Ri<r< R w 

n»|^> 

* H 

II 

o 


74 



Solid-fluid surfaces correspond to the heat pipe container internal wall 
interface with the working fluid, and the capillary wick interface with the fluid. 
The pipe wall conditions are no slip and matching temperature gradient. The 
conditions are: 


AXIAL RANGE 

RADIAL RANGE 

CONDITION 

x = 0 

R v < r < R/ 

u = 0 

v = 0 

dT l 

dx 

x = L 

R y <r<R l 

X ^ ^ * 

II II 

II O O 
o 

0 < x < L 

r-R, 

u = 0 
v = 0 

A- dTl -k dT " 
‘dr w dr 

0 <x< L 

r — R v 



The liquid-vapor interface of the heat pipe is the most difficult feature to 
analyze. The interface serves as the medium for communication between the 
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liquid and vapor along the entire axial length of the pipe. The flow charac- 
teristics are distinctly different between the respective liquid and vapor sides 
of the interface. The interface is subjected to physical phenomena through the 
action of surface tension that is generally not considered in the liquid or vapor 
flow fields . 2 The general liquid-vapor interface conditions are as the following: 
Kinematic Surface Condition 

Vh = 0 (60) 

where U lv is the velocity of the liquid-vapor interface surface. 

Conservation of Mass 

Pv^Ov-P/^OI (61) 

where V 0v is the vapor radial phase change velocity and F 0/ is the liquid radial 
phase change velocity. In future uses, V 0 will refer to the liquid phase change 
velocity unless for clarity V 0l is used. 

Continuity of Tangential Stress 

The tangential stress condition provides an important coupling relation be- 
tween the liquid and vapor phases. Since the vapor phase solution is not in- 
cluded in this study, as a first approximation the vapor tangential shear stress 
will be neglected and replaced with the axial velocity no slip condition 

2 The vapor and liquid regions are each assumed to be a single phase such 
that multi-phase phenomena in these bulk flow regions are not consid- 
ered. Entrainment of liquid by shearing action of the vapor, 
condensation in the bulk vapor, and vapor bubble formation in the liquid 
arc common examples of departures from the single phase assumption. 
These phenomena are outside the scope of the present analysis. 
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w 0 = ° 


(62) 


Continuity of Normal Stress 


P v ~ ?! + Pv ^ Ov - J/ 0 /) + 2 


ov[ ov v 

F v ~ — 

or or 


(63) 


Continuity of Thermal Flux 


, 37/ 

Pv 1 Ov h fg- - k l~T , — 
C r 


(64) 


where 




Z>, 


D p+Vw 


The term C w accounts for the effective liquid-vapor interface surface area due 
to the presence of the wick pores of diameter D p and of the screen wires of di- 
ameter D w . 

Supplemental Relations 

Supplemental relations can be used to derive relations between the liquid 
and vapor at the interface. Using the kinetic theory for interphase mass 
transfer [46] and the Clausius-Clapcvron Equation [47], the radial liquid ve- 
locity at the liquid-vapor interface is 


^0 = 
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(65) 
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The interphase heat flux can be expressed by using the continuity equation and 
substituting Eq. (65) in Eq. (64) 


*»- s/3r t s ) 


( 66 ) 


7.2 THE VAPOR PHASE MATHEMATICAL FORMULATION 

Heat pipe dynamics are strongly dependent on the vapor phase so that 
complete solution of the heat pipe problem requires solution of the full system 
of vapor phase governing equations. The intent of this research was to develop 
a computational algorithm for liquid phase dynamics that would be combined 
with a vapor phase algorithm. Development of the vapor phase algorithm 
encountered extreme difficulties such that a vapor code was not available for 
complete checkout of both the liquid and vapor codes at the time this work 
was done, see [44]. An approximation of the vapor phase was developed to 
permit demonstration of the liquid phase code operation. 

The characteristics of vapor flow allow some simplification in developing 
the approximate model. In general the vapor response is orders of magnitude 
faster than the liquid response such that the vapor can possibly be treated as 
quasisteady in a liquid time scale reference [31]. The vapor is also assumed 
to be uniform in a selected control volume. A mass balance on the control 
volume for a single arbitrary time step gives 



where M v is the total vapor mass at a given discrete time. S, is the time step 
size, and m v is the net mass flux due to evaporation and condensation. The 
sign of m Y is taken such that evaporation, which has a negative velocity, 
produces mass addition, while condensation produces mass extraction from the 
vapor space. 

Solution of Eq. (67) can be used to fix the vapor state. The solution is de- 
pendent on the phase change process, which is a function of the heat pipe op- 
erating limits. This study will consider restriction only due to the sonic limit. 
The solution is then dependent on whether the pipe operation is subsonic or 
sonic limited. 


7 . 2.1 Subsonic Vapor Operation 

Subsonic operation is assumed to occur such that the vapor is uniform at a 
given time step. The mass balance control volume is then the entire vapor 
space. The net phase change mass transfer can be expressed as 

m v= f* P v ^o v C xv 2t: R v cLx (68) 

J 0 

Using the continuity equation (Eq. (61)) and the kinetic energy formulation for 
V 0 (Eq. (64)) gives 


m x = C W R V 
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‘g 








(69) 
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where T v is the uniform vapor temperature, T s is the liquid-vapor interface 
surtacc temperature, and R s is the gas constant for the working fluid of inter- 
est. As a first approximation, all parameters except the temperatures in Eq. 
(69) are assumed to be at time level t< n l 

Substituting Eq. (69) in Eq. (67) and dividing by the fixed volume of the 
vapor space V gives 


P ( r" 



¥ 






(70) 


The Clausius-Clapevron Equation can be solved for the vapor density to give 


Pv = 



(71) 


All parameters except temperature are again assumed to be from time level 

t< n \ 

Equations (70) and (71) can be equated to produce a nonlinear expression 
for T y in terms of time level /<"> properties and time level r< n +>> liquid-vapor 
interface surface temperatures. The equation can be solved using Newton s 
method. 


SO 



7.2.2 Sonic Limit Operation 

The sonic limit condition is a restriction on axial mass and heat transfer. 
This restriction implies that the entire vapor space cannot be treated as uni- 
form. The vapor space will be approximated as two control volumes. The first 
control volume is the vapor space of the evaporator. The second control vol- 
ume consists of the combined vapor space of the adiabat and condenser re- 
gions. The two control volumes communicate through the axial vapor 
transport corresponding to sonic limited operation. The control volume mass 
balance for the evaporator is 

A/‘" + " = A/Jf - 6,m ve - 6, m sm (72) 

and for the adiabatic-condenser is 


M l? l) = M ic } ~ S i m vc + S t rh . son (73) 

where rh son is the sonic limit axial vapor flux. The heat flux associated with the 
sonic limit as a function of vapor properties at the point of sonic velocity is 
approximated by [8J 


Qson 


P V ^ 'Jg Ay C 

v 2(1 4- y ) 


(74) 


The sonic limit mass flux can be derived from the sonic limit heat transfer 
given by Eq. (74). Solving for mass transfer of sodium vapor produces 
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m son ^ \ , 
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(75) 
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The phase change mass flux for each control volume is found from Eq. (68) 
with limits of integration corresponding to the respective control volume. The 
evaporator control volume has phase change mass flux of 


m 


ve 


C L e 

= Pv^( 

Jfl 


Ov C w- * R v dx 


(76) 


or 


m ve ^ p yhjg 


L e T7 v '-T^- JV' 
J 0 


(77) 


Phase change mass flux for the adiabatic-condenser control volume is 


m v C = \ Pv V 0v C w 2 71 R v^ 


(78) 


or 



(L-L e )T7 v '-T7^ f L T/ix 

/ g 



(79) 


Equations (72) and (73) can be solved using Eqs. (75), (77) and (79) for the 
vapor density during sonic limit conditions for each of the control volumes. 
Using the Clausius-Clapeyron Equation as in the subsonic case produces a 
nonlinear expression for vapor temperature in each control volume that can 
be solved with Newton's method. 
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7.3 THE COMPUTATIONAL IMPLEMENTATION 


The governing equations, boundary conditions, property relations, and 
supplemental equations were coded for solution in a Fortran computer pro- 
gram. The calculations are performed at discrete grid points in the liquid and 
pipe wall. For clarity, the program First will be presented in a top level sense 
to establish the overall sequence of calculation. The detailed algorithm will 
then be discussed. Equations presented in other sections of the thesis are re- 
peated when discussed below for convenience. 

Overall organization of the program is shown in the flowchart in Figure 20. 
The algorithm begins by setting various parameters and initializing conditions. 
The time stepping section is then entered. Operating limits and property de- 
pendent parameters are updated. The input heat flux is updated as is the va- 
por regime. The iteration section is then entered. The first partial time step 
(f*) is solved for liquid and pipe temperatures, vapor temperature, and liquid 
velocities. The same steps are then repeated for the second partial time step 
(b"+!>). Iterations are repeated on the first and second time split equations 
until temperatures of the liquid surface and container wall have converged to 
a user specified tolerance. After time step convergence the thermophysical 
properties are updated. Heat transfer limits are determined. Input, output, 
evaporation, and condensation heat transfer rates are calculated. The 
capillary limit heat transfer is calculated and compared to evaporation heat 
transfer. A stop flag is set if the limit is exceeded. A steady state test is per- 
formed by comparing input and output heat fluxes. Flags are set if steady- 
state is reached. The liquid pressure gradient and liquid pressures are evalu- 
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atcd. The capillary radii are determined from vapor pressure, liquid pressure, 
and phase change velocities. Capillary radii are checked for the capillary 
pumping limit. A stop flag is set if the capillary limit is exceeded. Time step 
data are printed, and the program advances to the next time step unless a stop 
flag was issued or time has expired. If calculations are to be terminated, re- 
start data are printed to a restart file which can be used to continue the cal- 
culations in a later run. 


7.3. 1 Details of the Numerical Procedure 

The detailed program can be considered in three broad sections: the startup 
section, the iteration section, and the completion section. The startup section 
shown in Figure 21 begins the calculation cycle. Calculations can be per- 
formed beginning from a time zero condition or in restart mode to continue 
calculations of a previous case. The algorithm begins by setting common and 
dimension areas, and reading input data. Input data consist of the restart 
option flag, time step size, grid step sizes, vapor space radius, wick pore radius, 
screen wire radius, environment heat sink radiation temperature, and initial 
heat pipe temperature. The restart option flag controls the program initial- 
ization steps so that for a restart case, the program will initialize with the ap- 
propriate data from the previous case. Geometry parameters are calculated. 
Time parameters controlling the duration of the calculations are set. Data 
output print control parameters are set. Heat pipe surface radiation charac- 
teristics are set as are input heat flux control parameters used in Eqs. (SO) and 



(81). Thermophysical properties of the pipe wall material and sodium working 
fluid properties are initialized. Combinations of various constants used re- 
peatedly throughout the code are calculated to reduce total arithmetic. 
Peaceman-Rachford ADI parameters, operating limits and temperature and 
velocity fields are initialized. Startup data are printed, and control passes to 
the beginning of the time steps in the iteration section. 

The iteration section shown in Figure 22 begins the time step and performs 
iterations between temperatures and velocities until acceptable convergence is 
achieved. Property dependent operating limits and ADI parameters are up- 
dated. Input heat flux is updated for the new time step. The input heat flux 
profile and time dependence are user defined functions 

t < T 

given by Eqs. (80) and (81) 



and for 

t > T 

rqno(8 1 )q jn = <? max 

An option is included to smooth the abrupt step change of heat input at the 
evaporator-adiabatic boundary. The vapor operating regime is checked using 
Eq. (75). Operation is subsonic if axial heat transfer is below the sonic limit 
and is sonic limited otherwise. 
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Iterations are set up by saving data from the just completed time step and 
using these data as the initial guess for the next time step. Iterations are per- 
formed by first solving for time t * liquid and container temperatures. The fi- 
nite difference Peaccman-Rachford Alternating Direction Implicit method is 
used to solve the energy equation in two partial time steps. The method es- 
sentially involves solution by integrating the energy equation for a partial time 
.step (time t*) in the radial direction, and then integrating for the remaining 
time step (time r <n+l) ) in the axial direction. 

The first time split liquid temperatures are then used to find time t* vapor 
temperatures. Depending on whether the vapor is subsonic or sonic limited, 
control passes to the appropriate subroutine to use Newton's method to solve 
the vapor temperature equation. Subsonic vapor temperature is determined 
from Eqs. (70) and (71). Sonic limited vapor temperatures are determined 
from Eqs. (71), (72), and (73). Based on vapor temperatures and liquid-vapor 
interface surface temperatures the liquid phase change velocity is calculated 
from Eq. (65). The bulk axial velocity distribution is calculated from the con- 
tinuity equation and the liquid phase change velocity using the following 
equation. 



> 0 ( 1 ) + 2V v 0 (n) + V^i) 


(S3) 


Bulk axial velocity is transformed to a velocity profile using the following 
equation 
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u* = a + b r* + c r*~ + d r*' 


(S3) 


Finally, the radial velocity distribution is found from an equation given below 


v(A%r) = 
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Then, calculations for the First time split are completed. 

The second time split is set up by setting the temperatures and velocities for 
the First iteration of the (r <,!+1) ) time split equal to the time t* values. The liquid 
and container temperatures, vapor temperatures, and velocities for the second 
partial time step (r< ,J+l) ) are found as in the First time split with liquid and 
container temperatures found using the ADI equation 


■ 4471 - 


rt+1) 


+ 




ilftj-Btfj-l + Btfj+B, 




(85) 


The energy equation and the flow dynamics equations are tightly coupled 
at the liquid-vapor interface such that the solutions of each are interdependent. 
In addition, the equations contain nonlinear terms that required 
quasilinearization treatment, which requires iterations. Iterations are per- 
formed between temperatures and velocities until an acceptable convergence 
is achieved. 

The completion section of the program shown in Figure 23 includes all re- 
maining functions to complete the time step calculations and send control to 
the next time step or terminate execution. Thermophysical properties are up- 
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dated based on the temperature conditions. A mass balance is performed to 
track working fluid inventory. The sonic limit heat transfer is determined from 
Eq. (74). Capillary limit heat transfer is determined from 


A P 

Ax 


cap 


= - F iQ, 


cap 
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Input, output, evaporation, and condensation heat transfer rates arc evalu- 
ated. A stop flag is issued if the capillary limit is exceeded. Input heat transfer 
is compared with output heat transfer, and a steady state flag is issued if they 
arc approximately equal. 

Liquid pressures are found from the following equation which is found from 
the axial momentum integral equation 
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The integration is performed by assuming liquid pressure and vapor pressure 
are equal at the end of the condenser (.r = L ). The liquid-vapor interface ra- 
dial momentum condition given by Eq. (63) can be solved for the capillary 
radii of curvature to give 
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Figure 23: Algorithm Completion Flowchart 


The total liquid pressure drop is compared to the maximum possible capillary 
pumping capability of the wick. is given by 

AP ‘°p=jt (89) 
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A stop flag is activated if the limit is exceeded. Data from the time step are 
printed, and the next time step is initiated unless a stop flag has been activated 
or time for the case has expired. Data from the final time step arc dumped to 
a restart file, and execution is terminated. 


7.4 SOME RESULTS 

Capabilities of the computational algorithm were tested by running two 
heat pipe startup cases. The heat pipe dimensions were chosen to be consistent 
with the experimental device described by Merrigan et al. up to the limit of the 
uniform grid used in this analysis [27.28]. The experimental heat pipe in that 
study used lithium as the working fluid and an annular wick configuration. 
The annular region for liquid flow was formed between the pipe interior wall 
and a porous concentric tube constructed of 7.25 layers of pressed screen. The 
high temperatures involved required the container to be constructed from the 
refractory alloy molybdenum. The heat pipe geometry parameters used in this 
analytical study compared to the actual parameters are (see Figs. 18 and 19): 


Parameter 

Actual 

Approximate 

Heat Pipe Internal Length (L) 

4.0 

4.0 m 

Evaporator Length ( L e ) 

0.4 m 

0.4 m 

Condenser Length ( L c ) 

3.0 m 

3.0 m 

External Diameter (D H ) 

1.90 cm 

1.886 cm 

Internal Diameter {D t ) 

1.60 cm 

1.598 cm 

Vapor Space Diameter (D v ) 

1.49 cm 

' 1 .490 cm 
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Effective Wick Pore Diameter (D p ) 53 /im 


53 um 


where the effective wick pore diameter was experimentally determined from 
surface tension measurements. 

An axial grid step size of h x = 0.2 m and a radial grid step size of 
h r = 1.1 x 10~ 4 m were used. The computational grid consisted of 21 axial grid 
steps and 20 radial grid steps. The grid mapping to the heat pipe configuration 
is: 


Radial Liquid Space 

1 <j< 6 

Radial Pipe Wall 

6 < y < 20 


Axial 

Evaporator Region 

1 < / < 3 

Axial 

Adiabatic Region 

3 < / < 6 

Axial 

Condenser Region 

6 < i < 21 


Numerical experimentation showed that the code is highly sensitive to time 
step size. Too large a time step produced an inconsistent temperature at the 
liquid-vapor interface. The temperature was inconsistent in that the temper- 
ature gradient at the surface changed sign relative to the gradient of the inte- 
rior liquid. For example, in the evaporator energy is conducted from the pipe 
wall to the liquid-vapor interface. The liquid-vapor surface temperature 
should be the lowest temperature in the radial conduction path in order for the 
energy to reach the surface. Too large a time step will produce a physically 
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questionable surface temperature that is greater than the temperature of the 
interior liquid. Too large a time step was also found to cause oscillations in the 
pressure calculations. 

The time step size was selected by running a series of calculations for three 
hundred time steps each with gradually decreasing step sizes. A time step size 
of S,= 2 x 10 -3 second was found to eliminate inconsistent liquid-vapor inter- 
face temperatures, but did not completely remove pressure oscillations. A time 
step size of <5, = 1 x 10 -3 second removed the pressure oscillations for the three 
hundred time steps calculation. This time step was used in the calculations. 


7 . 4.1 Slow Startup Test 

A startup test was devised to supply input energy to the evaporator such 
that the heat pipe would warm up without reaching the sonic limit. This test 
was to check overall program implementation. Input heat flux was supplied 
according to the function 
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t < T 



t > T 


9 in ^max (91) 

where <7 max is the maximum input heat flux and t is the time stretching pa- 
rameter that controls how quickly q in reaches q max . In this test. 

*7max =1x10° IVjm 2 and z - 300 seconds corresponding to a five minute pe- 
riod for the input heat flux to reach the maximum. 

Calculations were performed for 1 14,193 time steps, or a period slightly 
greater than 1 14 seconds. The run was terminated at this time since the liquid 
pressure drop exceeded the capillary pumping limit. The liquid pressure drop 
and the capillary pumping limit as functions of time are shown in Figure 24. 
Analysis of the results shows that the capillary pumping limit was exceeded 
due to a large jump in pressure drop at a single time step. Figure 24 also 
shows that the pressure calculations had been experiencing oscillations with 
usually small amplitudes of approximately 2.5 Pa. The oscillations began after 
the 300 time steps used to verify the adequacy of the time step size of 
1 x 10 -3 second. The liquid pressure drop curve is drawn through the centers 
of the oscillation ranges in Figure 24. The oscillations become very large at 
t = 90 seconds, but begin to decay until the large jump occurs at approximately 
1 1 4 seconds. 


96 



1 0 


CO 

CL 


CL 

O 


<D 


D 

CO 

CO 

0 ) 



Figure 24: Slow Startup Pressure Dr 
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in detail in Figure 25. The liquid 
i consecutive time steps for 



19,990 < k < 20,000, corresponding to the time around time t = 20 seconds in 
Figure 24. The oscillations arc of essentially constant amplitude and well be- 


haved until time t > 80 seconds as shown in Figure 24. 
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Figure 25: Slow Startup Liquid Pressure Drop Oscillations 


Power transfer as a function of time is shown in Figure 26. While the input 
power briefly exceeds the sonic limit line, the phase change heat transfer re- 
mains below the sonic limit. At approximately 103 seconds into the transient, 
the sonic limit exceeds the capillary limit so that the capillary limit becomes the 
controlling heat transport limit. At approximately 105 seconds, the unlikely 
result is produced that the phase change heat transfer exceeds the input heat 
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transfer. With the exception of this anomaly, the heat transfer curves arc well 
behaved, smooth functions of time. In addition, the phase change heat trans- 
fer is well below the capillary limit heat transfer at the point when liquid 
pressure drop jumps to exceed the capillary pumping limit. The cause of the 
liquid pressure drop jump is not apparent from the heat transfer curves. 

The vapor temperature as a function of time is shown in Figure 27. The 
vapor temperature is also a smooth, although rapidly increasing, function of 
time. The cause of the jump in liquid pressure drop is also not apparent from 
vapor temperature. 

The pressure drop calculation can be seen in Figure 24 to produce 
oscillatory behavior with large sporadic oscillations for time greater than 80 
seconds. The pressure calculations are a sensitive function of time step size. 
Numerical experimentation also showed that the pressure calculations arc 
sensitive to property updates and changes in the liquid-vapor interface 
boundary conditions corresponding to the change between subsonic and sonic 
limited operation. Sensitivity to property updates was identified bv varying 
the frequency at which updated properties are calculated. 

There are many potential causes for the pressure oscillations. The pressure 
oscillations may be physical, but this possibility is remote considering the nu- 
merical sensitivities of the method. The time step size of I x 10~ 3 second may 
have been too large. If this is the case, a smaller time step may not be ac- 
ceptable considering the cost of computations for a transient process that oc- 
curs over a period of tens of minutes. The oscillations may have resulted from 
the derivation of the axial velocity profile in which the shape is determined by 
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Figure 26: Slow Startup Heat 
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uid pressure occurs at the end of the condenser, although this assumption is 
not generally valid for the wide range of operating conditions encountered 
during a startup transient. 





7.4.2 Sonic Limit Startup Test 

A startup test was performed to test algorithm calculation of sonic limited 
operation. Using again r = 300 seconds, the maximum input heat flux was 
increased to q miX = 1.5 x 10 6 Wj m 2 to assure the phase change heat transfer 
would reach the sonic limit. Heat transfer results are shown in Figure 28. 

Input power exceeds the sonic limit at 20 seconds, and remains above the 
sonic limit through the duration of the test of 130 seconds. The evaporation 
and condensation heat fluxes are essentially equal as the sonic limit is reached 
at approximately 40 seconds. Evaporation and condensation follow the sonic 
limit curve until 80 seconds. After 80 seconds, condensation heat flux in- 
creases above the sonic limit, while evaporation heat flux continues to follow 
the sonic limit curve. This unlikely result is a result of the inadequacies of the 
simple vapor model used to perform calculations. The vapor model uninten- 
tionally forces a high condensation rate, which forces the condenser temper- 
ature to increase rapidly due to the high thermal resistance of the radiation 
boundary condition on the condenser. This temperature effect is shown in 
Figure 29. The evaporator vapor temperature smoothly and gradually in- 
creases with time, while the condenser vapor temperature rapidly increases. 
The condenser temperature eventually increases sufficiently so that radiation 
heat transfer from the condenser becomes very efficient. At approximately 93 
seconds, the output heat transfer exceeds the evaporation heat transfer, which 
is another unlikely result. Calculations are terminated at 130 seconds due to 
the condensation heat transfer exceeding the capillary heat transfer limit. 

Vv hile the capillary limit is not strictly a limit on condensation, calculations 
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7.5 CONCLUSIONS AND RECOMMENDATIONS 

The heat pipe is a very complex device such that performing analysis of heat 
pipe transients is a difficult venture. The difficulty is due to the coupled liquid 
and vapor dynamics, the phase change process, the coupled heat transfer 
problem with nonlinear boundary conditions, the physical geometry, and the 
transient time requirement. This study has presented the full system of gov- 
erning equations, including boundary conditions, required to solve the liquid 


104 




phase heat pipe problem. A simplified solution was formulated by using cer- 
tain assumptions and integrating the liquid dynamics equations. 

The resulting formulation identified a key requirement for future research 
in solving the heat pipe problem. This study used a kinetic theory approach 
to model the phase change process as expressed by Eq. (65). The model con- 
sists of a large coefficient multiplying a very small temperature difference. The 
model requires a very small time step for stability even with the simplifications 
used in this analysis. The large coefficient also transforms temperature differ- 
ences that are otherwise beyond machine accuracy into computationally sig- 
nificant terms. Future research should investigate an alternative model for the 
phase change process. 

Difficulty with geometry was encountered due to the widely disparate 
length scales in the different coordinate directions. The radial dimension 
across the liquid is much smaller than the radial dimension across the pipe 
wall, which in turn is much smaller than dimensions in the axial direction. 
One approach to address this problem is to assume radial gradients are negli- 
gible. The radial pressure gradient was neglected in this study. This study also 
found small radial temperature gradients indicating that the radial thermal 
resistance is small. The drawback of this approach is the loss of fidelity with 
the true physics of the problem. A preferable approach may be to use a vari- 
able mesh grid with a to-be-determined simplified treatment of radial gradi- 
ents. The variable mesh grid allows fine resolution in important regions with 
coarse resolution in the other regions. The simplified treatment of radial gra- 
dients would take advantage of the small radial gradients. 



The long duration of transients presents a competing concern with model 
complexity. An overly complex model with high fidelity may be too computa- 
tionally expensive for practical calculations of transients. 

This study has also shown that the transient heat pipe problem is not ame- 
nable to a straightforward simplification such as is used for the vapor phase. 
The solution of the heat pipe transient problem requires a full solution of the 
liquid and vapor phases. This elusive solution is left to future efforts. 
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